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Prof.  Geoffrey  M.  Lilley  to  Join  ICASE 

Professor  Geoffrey  M.  Lilley,  a  world  renown 
fluid  dynamicist  and  an  expert  on  turbulence 
and  acoustics,  will  join  ICASE  in  April  1998  as 
Chief  Scientist  for  Fluid  Mechanics.  In  this  posi¬ 
tion,  Prof.  Lilley  will  have  considerable  indepen¬ 
dence  to  expand  ICASE’s  aeroacoustics  program 
through  interactions,  collaborations,  and  mentor¬ 
ship  of  ICASE  staff  scientists,  ICASE  consultants, 
as  well  as  NASA  and  industry  research  scientists. 

Prof.  Lilley  obtained  the  B.Sc.  degree  in  En¬ 
gineering  with  First  Class  Honors  from  Imperial 
College  in  1943.  Later,  in  1945,  under  the  super¬ 
vision  of  Professors  Sir  Leonard  Bairstow  and  Sir 
George  Temple,  Prof.  Lilley  obtained  the  M.  Sc. 
and  D.I.C.  degrees  also  from  Imperial  College. 

In  October  1946,  the  College  of  Aeronautics  at 
Cranfield  was  established  in  an  airfield  used  by  the 
RAF.  The  training  at  Cranfield  was  to  be  at  the 
post  graduate  level  with  research  to  be  undertaken 
by  staff  and  students.  Prof.  Lilley  joined  the  Col¬ 
lege  as  a  founding  member  and  lectured  until  1955. 
In  1955,  he  was  appointed  Deputy  Head  of  the  Col¬ 
lege  of  Aeronautics.  Prof.  Lilley  held  this  position 
for  8  years  until  he  joined  the  University  of  South¬ 
hampton  as  Professor  and  Head  of  the  Aeronau¬ 
tics  and  Astronautics  Department  in  1963.  He  re¬ 
mained  Department  Head  for  15  years.  In  1983, 
Prof.  Lilley  became  Professor  Emeritus. 

Prof.  Lilley’s  long  career  has  provided  him  with 
a  wealth  of  knowledge  and  experience  in  both  the¬ 
oretical  and  experimental  aspects  of  fluid  mechan¬ 
ics.  He  has  been  responsible  for  the  design,  man¬ 
ufacture,  and  installation  of  many  wind  tunnels  at 
Cranfield.  He  has  also  published  and  made  im¬ 


portant  contributions  in  the  fields  of  aeroacoustics, 
flutter,  sonic  boom,  turbulence,  aerodynamics  of 
road  vehicles,  human-powered  flight,  theory  of  jets, 
propeller-wing  interference,  and  the  design  of  su¬ 
personic  civil  transport  and  AVSTOL  aircraft. 

Among  the  many  honors  he  has  received,  partic¬ 
ularly  worth  mentioning  are  the  Gold  Medal  of  the 
Royal  Aeronautical  Society  in  1983  and  the  AIAA 
Aeroacoustic  Medal  in  1984. 


Professor  Geoffrey  Lilley 
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ICASE:  In  Retrospect 


Amiram  “Ami”  Harten  David  Kamowitz  Milton  E.  “Milt”  Rose 

9/22/46-8/5/94  10/15/58-9/16/89  5/22/25-8/22/93 

Visiting  Scientist,  Staff  Scientist  (’86-’89)  Director  (’77-’86) 

Consultant  (76- ’94) 

Message  from  the  Director,  Manuel  D.  Salas 

The  success  ICASE  has  had  over  the  last  25  years  in  advancing  research  in  applied  mathematics,  computer 
science,  and  fluid  mechanics  has  been  a  result  of  the  vision  of  its  former  directors,  the  excellence  of  its 
scientific  staff  and  consultants,  and  the  dedication  and  commitment  of  its  administrative  staff.  We  dedicate 
this  issue  of  the  quarterly  to  the  memory  of  Ami  Harten,  David  Kamowitz,  and  Milt  Rose.  They  are  greatly 
missed  by  their  friends  and  colleagues;  their  contributions  to  science  live  on. 

Ami  Harten  was  a  brilliant  scientist.  He  pursued  his  doctoral  studies  under  Peter  Lax  at  the  Courant 
Institute  of  Mathematical  Sciences.  His  short  career  was  extremely  influential  to  those  involved  in  the 
numerical  solution  of  hyperbolic  conservation  laws.  Among  his  many  contributions  are  his  work  on  Total- 
Variation-Diminishing  schemes  and  Essentially  Non-Oscillatory  schemes.  Ami  was  a  cheerful,  outgoing 
fellow,  always  willing  to  share  his  insights  and  engage  in  friendly  discussion.  His  visits  to  ICASE  for  a 
period  covering  18  years  always  enriched  the  atmosphere  of  the  institute. 

David  Kamowitz  came  to  ICASE  in  1986.  His  Ph.D.  work,  under  Seymour  Parter  at  the  University 
of  Wisconsin,  was  on  multigrid  methods  for  elliptic  boundary  value  problems.  He  was  diagnosed  with  a 
rare  form  of  cancer  in  the  spring  of  1987.  Throughout  his  two  and  a  half  year  ordeal,  he  was  incredibly 
optimistic  and  courageous.  Despite  his  illness,  he  made  contributions  to  the  numerical  implementation  of 
outflow  boundary  conditions  and  the  application  of  multigrid  to  singular  perturbation  problems. 

Milton  E.  Rose  had  probably  the  greatest  influence  in  setting  ICASE’s  future  course.  I  met  him  when 
he  became  ICASE  Director  in  1977.  Milt,  as  he  liked  to  be  called,  was  a  career  administrator  with  a  love 
for  mathematics  and  a  sharp  sense  for  identifying  young,  talented  scientists.  His  doctoral  advisor,  Richard 
Courant,  taught  him  the  importance  of  understanding  and  incorporating  physical  ideas  into  computational 
mathematics.  His  interest  in  computational  fluid  dynamics,  his  gentle,  modest  nature,  and  his  disdain  for 
bureaucracy  was  an  instant  hit  with  me  and  we  remained  close  friends  until  his  death  in  1993.  Hanging 
on  my  wall  is  a  letter  he  wrote  to  me  on  April  8,  1981.  His  enthusiasm  for  carrying  out  research  and  his 
humble  nature  rings  through  it.  His  signature  now  fading,  it  reads: 

Dear  Manny: 

Our  experiments  have  been  sidetracked. . .  Simply  stated,  these  arose  from  an  unsuccessful  attempt 
to  impose  dissipative  boundary  conditions  at  all  boundaries. . . 

...  I  am  focusing  upon  treating  some  simpler  problems. . .  Hopefully,  with  these  clarified,  a  more 
focused  attempt  to  understand  the  remaining  problems  can  be  undertaken. 

Friedrichs  once  complained  to  me  that  applied  mathematicians  (he  included)  seem  to  engage 
themselves  in  clarifying  methods  which  scientists  and  engineers  have  already  established!  I’m 
reminded  of  the  story  of  the  fellow  who  enjoyed  show  business  -  his  job  was  sweeping  the  excrement 
[Milt  used  a  shorter,  more  descriptive  word]  of  the  elephants  in  the  circus!  Sometimes  I  think  of 
this  when  I  attempt  to  catch  up  to  where. .  .others  are  (and  have  already  passed!). 

Well,  we  also  serve  who  only  sweep! 

Milt. 
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ICASE  Administrative  Staff 


Front  row,  1-to-r:  Gwendolyn  W.  Wesson  -  Contract  Accounting  Clerk  (1994);  Linda  T.  Johnson  -  Office 
and  Financial  Administrator  (1974) 

Back  row:  Barbara  A.  Cardasis  -  Administrative  Secretary  (1986);  R.  Anne  Lomas  -  Payroll  and  Ac¬ 
counting  Clerk  (1992);  Etta  M.  Blair  -  Accounting  Supervisor  (1984);  Emily  N.  Todd  -  Conference  Manager 
(1982);  Shelly  M.  Johnson  -  Executive  Secretary /Visitor  Coordinator  (1989) 


Current  ICASE  Staff 

(ANM  =  Applied  and  Numerical  Mathematics;  CS  =  Computer  Science;  FM  =  Fluid  Mechanics) 


Research  Fellows: 

Staff  Scientists: 

Mavriplis,  Dimitri 

(ANM) 

Allan,  Brian  G. 

(ANM) 

Mehrotra,  Piyush 

(CS) 

Arian,  Eyal 

(ANM) 

Guattery,  Stephen  M. 

(CS) 

Associate  Research  Fellow: 

Interrante,  Victoria  L. 

(CS) 

Keyes,  David  E. 

(CS  &  ANM) 

Luo,  Li-Slii 

(FM) 

Ma,  Kwan-Liu 

(CS) 

Senior  Staff  Scientists: 

Povitsky,  Alexander 

(FM) 

Crockett,  Thomas  W. 

(CS) 

Girimaji,  Sharath  S. 

(FM) 

System  Administrators: 

Lewis,  R.  Michael 

(ANM) 

Clancy,  Leon  M. 

Loncaric,  Josip 

(ANM) 

Hess,  Bryan  K. 

Ristorcelli,  J.  Ray 

(FM) 

Sidilkover,  David 

(ANM) 

Zhou,  Ye 

(FM) 

ISE  currently  has  14  graduate  students  and  three  visiting  scientists. 

ICASE  Directors 

Facts  About  ICASE 

James  M.  Ortega 

1972-1977 

ICASE  Technical  Reports  to 

Date 

1,440 

Milton  E.  Rose 

1977-1986 

Conferences  /Workshops /Short  Courses 

55 

Robert  G.  Voigt 

1986-1991 

Hardbound  Volumes  Published  (since  1992) 

20 

M.  Yousuff  Hussaini 

1992-1996 

(Acting  Director  1991-1992) 

Manuel  D.  Salas 

1996-present 
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3D  Flow  Visualization  Using  Volume  Line 
Integral  Convolution 

Victoria  Interrante 1  and  Chester  Grosch 2 

Introduction 

Line  integral  convolution  (LIC)  is  a  flow-driven 
texture  generation  method  that  has  become  one 
of  the  best-known  and  most  commonly  used  tech¬ 
niques  in  computer  graphics  for  visualizing  2D  flow, 
or  flow  over  a  surface  in  3D.  The  popularity  of  LIC 
as  a  tool  for  3D  flow  visualization,  or  the  depiction 
of  flow  through  a  volume,  has  been  relatively  lim¬ 
ited  in  contrast,  however,  primarily  due  to  the  dif¬ 
ficulties  inherent  in  clearly  and  effectively  portray¬ 
ing  a  dense  volume  texture  in  a  static,  2D  image. 
Over  the  past  months,  we  have  been  investigating 
strategies  for  more  effectively  using  3D  LIC  for  the 
visualization  of  3D  flow.  Much  of  this  work  is  de¬ 
scribed  in  our  ICASE  Report  No.  97-35.  In  this 
article  we  highlight  new  results  from  our  continuing 
work  in  this  area. 

Background  and  Motivation 

Given  a  vector  field  and  an  input  texture,  line 
integral  convolution  produces  an  output  texture 
in  which  the  data  values  are  highly  correlated  in 
the  direction  of  the  flow.  Our  work  focuses  on 
methods  for  effectively  representing  the  flow  infor¬ 
mation  contained  in  the  dense  volumetric  textures 
produced  by  3D  LIC.  Our  strategies  include  selec¬ 
tively  emphasizing  flow  information  in  critical  re¬ 
gions  of  interest  in  the  volume  and  clarifying  the  3D 
structure  of  the  flow  by  facilitating  the  perceptual 
differentiation  of  the  densely  clustered  streamlines. 

Region  of  Interest  Definition 

By  concentrating  the  3D  texture  in  the  most  sig¬ 
nificant  areas  of  the  flow,  we  can  clarify  the  visual 
representation  of  the  data  and  facilitate  the  appre¬ 
ciation  of  the  most  relevant  information. 

When  LIC  is  used  in  conjunction  with  a  region 
of  interest  (ROI)  definition  based  on  the  value  of  a 
scalar  quantity  across  the  volume,  we  have  found 
that  best  results  are  achieved  when  the  ROI  mask 
is  applied  to  the  input  texture  rather  than  to  the 
output.  Figure  1  compares  the  two  effects. 

xDr.  Interrante  received  her  Ph.D.  from  the  University  of 
North  Carolina  at  Chapel  Hill,  and  joined  ICASE  as  a  Staff 
Scientist  in  1996. 

2 Dr.  Grosch  is  a  Professor  of  Oceanography  and  Com¬ 
puter  Science  at  Old  Dominion  University.  He  has  been  a 
consultant  at  ICASE  since  1980. 


Figure  1:  Clockwise  from  upper  left:  a  2D  slice  from 
a  3D  LIC  texture;  a  region  of  interest  mask,  defined 
by  velocity  magnitude;  results  when  the  ROI  mask  is 
applied  to  the  texture  generated  by  LIC;  results  when 
the  ROI  mask  is  applied,  as  a  preprocess,  to  the  input 
texture  whose  values  are  then  convolved  by  LIC. 

When  the  ROI  mask  is  applied  as  a  post  process, 
the  visibility  of  the  flow  information  is  directly  de¬ 
fined  by  the  values  in  the  ROI,  whose  edges  are  not, 
in  general,  guaranteed  to  be  everywhere  tangent  to 
the  direction  of  the  flow.  When  the  ROI  mask  is 
applied  as  a  preprocess,  the  same  basic  segmenta- 


Figure  2:  Upper:  ridges  of  velocity  magnitude  can  be 
used  to  define  a  surface  of  interest  in  the  flow;  Lower: 
3D  LIC  applied  to  a  texture  of  points  evenly  distributed 
over  this  velocity  ridge  surface. 
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tion  is  in  effect,  however  the  flow  itself  is  allowed 
to  define  the  explicit  boundaries  of  the  visible  in¬ 
formation. 

Figure  2  illustrates  a  second  method  for  ROI  def¬ 
inition.  Here,  a  ridge-finding  algorithm  is  used  to 
define  a  precise  surface  of  interest  in  the  volume, 
and  the  LIC  texture  is  derived  from  a  set  of  Gaus¬ 
sian  spots  uniformly  distributed  over  this  surface. 
However,  the  directional  information  is  not  pro¬ 
jected  onto  the  2D  surface.  All  calculations  are 
done  in  3D,  so  that  the  tufts  in  the  output  texture 
will  accurately  reflect  the  local  3D  orientation  of 
the  flow  in  the  immediate  vicinity  of  the  specified 
surface  of  interest. 

Clarifying  the  Dense  Texture  Data 

When  the  LIC  output  densely  occupies  a  3D  re¬ 
gion  of  space,  individual  streamlines  can  be  diffi¬ 
cult  to  discriminate  due  to  their  similar  shading. 
The  three-dimensional  spatial  relationships  among 
the  overlapping  streamlines  represented  in  the  LIC 
texture  volume  will  be  clarified  if  the  depth  discon¬ 
tinuities  in  the  scene  are  emphasized  through  the 
use  of  ‘visibility-impeding  halos’,  as  demonstrated 
in  Figure  3. 
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Figure  3:  Upper  left:  volume  rendered  streamlines  in 
a  flow  through  a  circular  jet  with  tabs;  Upper  right: 
visibility-impeding  halos  emphasize  depth  discontinu¬ 
ities;  Bottom:  haloed  streamlines  in  a  flow  through  a 
rectangular  aperture. 

Information  about  the  forward/backward  direc¬ 
tion  of  the  flow  can  be  incorporated  into  a  static 
2D  image  through  the  use  of  tapered  streamlines. 


Figure  4:  Tapered  streamlines  convey  the  for¬ 

ward/backward  direction  of  the  flow. 


A  simple  modification  to  the  fast- LIC  algorithm  al¬ 
lows  the  efficient  computation  of  oriented  stream¬ 
lines,  as  shown  in  Figure  4.  This  modification  is 
based  on  the  use  of  an  asymmetric  filter  of  the  form: 

/ length 

I'd  -  Y  c  < 1 

i= 0 

in  place  of  a  simple  box  filter.  Texture  values  at 
subsequent  voxels  along  a  streamline  may  still  be 
computed  incrementally  based  on  previous  values: 

I\  —  {I’o  -  vq)/c  +  V f  length +  \clle"9th 

I-l=(Io-VfUngthCflen9th)c+V 

Future  Work 

We  are  continuing  to  investigate  methods  for 
more  efficiently  generating  smooth,  cyclic  anima¬ 
tions  of  3D  LIC  textures  along  streamlines  in 
steady  flow  data,  and  are  actively  working  on  ex¬ 
tending  our  3D  LIC  algorithm  to  the  visualization 
of  streaklines  in  3D  unsteady  flow  data. 

References 

V.  Interrante  and  C.  Grosch,  “Strategies  for  Ef¬ 
fectively  Visualizing  3D  Flow  with  Volume 
LIC,”  Proceedings  of  IEEE  Visualization  ’ 97 ' 
pp.  421  -  424. 

C.  E.  Grosch,  J.  M.  Seiner,  M.  Y.  Hussaini  and 
T.  L.  Jackson,  “Numerical  Simulation  of  Mix¬ 
ing  Enhancement  in  a  Hot  Supersonic  Jet,” 
Physics  of  Fluids,  9(4):  1125  -  1143,  1997. 

V.  Interrante,  “Illustrating  Surface  Shape  in 
Volume  Data  via  Principal  Direction- Driven 
3D  Line  Integral  Convolution,”  Computer 
Graphics  Proceedings ,  Annual  Conference 
Series ,  pp.  109  -  116. 


5 


Eigenvalue  Bounds  Using  Graph 
Embeddings 

Stephen  Guattery 1 

Introduction  and  Motivation 

Combinatorial  techniques  provide  a  useful  way 
of  bounding  extreme  eigenvalues  of  Laplacian  and 
related  matrices.  Any  symmetric  matrix  can  be 
described  in  terms  of  a  weighted  underlying  graph. 
To  construct  bounds,  one  first  constructs  embed¬ 
dings  of  pairs  of  these  graphs  into  each  other,  then 
derives  quantities  that  characterize  the  embedding. 
Bounds  are  computed  from  these  quantities.  Em¬ 
bedding  methods  have  the  advantage  that  they  can 
often  be  applied  to  classes  of  matrices:  specify¬ 
ing  the  form  of  the  embedding  for  the  class  yields 
bounds  for  all  members  of  the  class.  While  these 
techniques  were  originally  developed  for  Laplacian 
matrices,  a  number  of  recent  developments  have  ex¬ 
tended  them  to  more  general  classes  of  symmetric, 
positive  definite  matrices.  We  review  some  bound¬ 
ing  techniques  and  discuss  how  they  have  been  ex¬ 
tended. 

Bounds  on  the  extreme  eigenvalues  have  many 
uses.  The  Laplacian  is  used  in  representing  phys¬ 
ical  problems  involving  elliptic  partial  differential 
equations.  Since  these  matrices  are  symmetric, 
their  extreme  eigenvalues  can  be  used  in  computing 
their  spectral  condition  numbers,  which  are  used  to 
bound  the  rates  of  convergence  of  iterative  linear 
system  solvers,  and  to  analyze  the  quality  of  pre¬ 
conditioners. 

The  smallest  nonzero  eigenvalue  A2  of  a  Lapla¬ 
cian  is  related  to  properties  of  the  associated  graph. 
This  connection  has  been  exploited  in  the  design 
and  analysis  of  graph  algorithms,  particularly  al¬ 
gorithms  for  finding  small  separators:  Bounds  on 
A2  can  be  used  to  compute  bounds  on  cut  quality 
and  to  isolate  the  structure  of  the  eigenvectors  on 
which  the  cuts  are  based.  The  eigenvalue  A2  is  also 
related  to  expansion  properties  of  graphs,  and  can 
be  used  to  determine  if  a  graph  is  an  expander. 

Laplacians,  Graphs,  and  Embeddings 

A  Laplacian  matrix  can  be  defined  in  terms  of 
a  graph  as  follows:  Let  G  be  an  undirected  graph 
with  vertices  ui, . . .  ,un.  Then  the  Laplacian  of  G 

1Dr.  Guattery  received  his  Ph.D.  from  Carnegie  Mellon 
University.  He  joined  ICASE  as  a  Staff  Scientist  in  1995. 


is  an  n  x  n  matrix  L  such  that 

degree  (u^)  if  i  —  j 
—  1  if  (i,  j)  is  an  edge  of  G 

0  otherwise. 

Laplacians  are  symmetric  and  have  all  row  sums 
equal  to  zero.  We  assume  for  the  discussion  below 
that  the  matrix  is  irreducible,  which  corresponds 
to  assuming  that  G  is  connected. 

The  definition  can  be  extended  to  weighted 
graphs  if  the  nonzero  off-diagonal  entries  have  neg¬ 
ative  values  other  than  —1;  these  entries  are  the 
negatives  of  the  weights  of  the  edges  in  the  graph. 
In  this  case  we  define  the  degree  of  a  vertex  to  be 
the  sum  of  the  weights  of  its  incident  edges.  In 
either  the  unweighted  or  weighted  case,  if  the  row 
sums  are  not  equal  to  zero,  we  can  think  of  the 
surplus  (or  deficiency)  of  the  sum  for  row  i  as  rep¬ 
resenting  the  weight  of  an  edge  to  a  zero  Dirich- 
let  boundary,  which  can  be  represented  as  a  single 
ground  vertex  in  the  graph.  This  vertex  is  implicit 
in  the  matrix  structure. 

The  techniques  discussed  below  make  use  of 
graph  embeddings:  Let  G  and  H  be  connected 
graphs  on  the  same  set  of  vertices.  An  embedding 
of  G  into  H  is  a  set  of  paths  in  H  such  that  for 
every  edge  (i,  j)  in  G,  the  set  includes  a  path  from 
vertex  i  to  vertex  j. 

Note  that  an  edge  in  H  could  end  up  on  more 
than  one  path.  To  capture  this  possibility,  we  in¬ 
troduce  the  idea  of  congestion .  Congestion  is  sim¬ 
ple  in  the  unweighted  case:  the  congestion  of  an 
edge  e  in  H  is  the  number  of  paths  in  the  embed¬ 
ding  that  include  e.  In  the  weighted  case,  we  assign 
each  path  the  weight  of  the  corresponding  edge  in 
G.  The  congestion  of  an  edge  e  in  H  is  the  sum 
of  the  weights  of  the  paths  that  include  e  divided 
by  the  weight  of  e.  The  congestion  of  a  path  is  the 
sum  of  the  congestions  of  the  edges  in  the  path. 

The  Path  Resistance  Method 

Using  the  ideas  of  graph  embeddings  and  path 
resistances,  Guattery,  Leighton,  and  Miller  have 
shown  a  way  to  compute  a  lower  bound  on  the 
smallest  nonzero  eigenvalue  A2  of  a  Laplacian.  Let 
G  be  the  underlying  graph  of  the  Laplacian  L\  for 
simplicity,  we  assume  that  G  is  unweighted.  Let 
K  be  the  complete  graph  on  the  n  vertices  of  G  (a 
complete  graph  has  an  edge  between  every  pair  of 
vertices).  The  bound  is  based  on  an  embedding  of 
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K  into  G .  The  technique  is  called  the  path  resis¬ 
tance  method ,  and  is  described  as  follows: 

1.  Construct  a  path  embedding  of  K  into  G. 

2.  Compute  the  congestion  Cjj  of  each  edge  (i,j) 
in  G. 

3.  For  each  path  P,  allocate  a  resistor  of  size  Cij 
to  edge  Cij  on  P. 

4.  For  each  path  P  compute  its  resistance  as  a 
series  resistive  circuit,  i.e.,  Y^eij^Pcij’  Let  r 
be  the  maximum  resistance  over  all  paths. 

5.  return  “n/r  <  A2”. 

The  proof  that  this  procedure  gives  a  lower 
bound  is  based  on  decomposing  the  problem  into 
edges  from  K  and  corresponding  paths  in  G .  Each 
edge  in  each  path  in  G  is  scaled  by  its  congestion. 
The  resulting  path  problems  have  an  interpretation 
in  terms  of  the  theory  of  resistive  circuits,  hence  the 
name  “path  resistance  method.” 

The  path  resistance  method  can  be  extended  in 
a  number  of  ways.  Some  of  the  extensions  allow 
its  application  to  a  wider  range  of  graphs  including 
weighted  graphs  and  diagonally  dominant  graphs 
with  positive  off-diagonal  entries.  Another  exten¬ 
sion  lets  the  method  be  used  with  Laplacians  with 
zero  Dirichlet  boundary  conditions;  this  variant 
embeds  a  star  rather  than  a  complete  graph.  Fi¬ 
nally,  priorities  can  be  assigned  to  paths  in  the  em¬ 
bedding;  these  priorities  affect  the  relative  amounts 
of  congestion  produced  by  various  paths,  and  can 
be  used  in  certain  cases  to  give  better  lower  bounds. 

The  best  bounds  produced  by  the  path  resistance 
method  using  priorities  are  not  tight.  This  is  a 
consequence  of  a  result  by  Nabil  Kahale.  However, 
Guattery,  Leighton,  and  Miller  have  shown  that 
for  the  Laplacian  of  any  unweighted  tree  T,  the 
path  resistance  method  as  stated  above  (i.e.,  using 
uniform  priorities)  gives  a  lower  bound  that  is  off 
by  at  most  a  factor  of  0(logdiameter(T)). 
Embeddings,  Generalized  Eigenvalues,  and 
Spectral  Condition  Numbers 

In  his  Ph.D.  thesis,  Keith  Gremban  used  similar 
techniques  for  bounding  spectral  condition  num¬ 
bers  of  preconditioned  systems.  In  particular,  if 
A  and  B  are  positive  definite  Laplacians,  he  has 
shown  how  to  use  embedding  techniques  to  give  an 


upper  bound  on  k(B~1A),  where  k  is  the  spectral 
condition  number.  Gremban ’s  techniques  require 
embedding  B  into  A  and  A  into  B. 

Gremban  used  the  results  to  bound  the  condition 
number  of  Support  Tree  Preconditioners,  which 
are  tree-structured  preconditioners  built  on  an  ex¬ 
tended  vertex  set  (i.e.,  the  matrix  for  the  precondi¬ 
tioner  has  a  somewhat  larger  order  than  the  orig¬ 
inal  matrix).  Because  of  the  tree  structure,  these 
preconditioners  are  amenable  to  fast  parallel  com¬ 
putation  in  preconditioned  conjugate  gradient  cal¬ 
culations. 

Gremban  gives  extensions  to  cases  where  A  and 
B  may  have  positive  off-diagonal  entries,  but  re¬ 
quires  that  the  matrices  be  diagonally  dominant. 
Guattery  has  given  an  extension  to  non-diagonally 
dominant  positive  definite  matrices,  which  allows 
for  incomplete  factor  preconditioners. 

We  have  successfully  applied  these  techniques  to 
a  wide  variety  of  preconditioners  including  Sup¬ 
port  Tree  preconditioners,  preconditioners  based 
on  incomplete  Cholesky  factorization,  block  diag¬ 
onal  domain  decomposition  preconditioners,  and 
preconditioners  for  various  basic  iterative  methods. 

An  Exact  Relationship  Between  Embed¬ 
dings  and  Eigenvalues 

As  noted  above,  the  bounds  produced  by  these 
embedding  methods  are  not  tight.  However,  there 
is  a  fundamental  connection  between  embeddings 
and  graph  spectra.  Guattery  and  Miller  have 
shown  that,  by  changing  the  representation  of  the 
embedding  to  include  (arbitrary)  edge  directions, 
it  is  possible  to  use  a  specific  embedding  to  pro¬ 
duce  a  matrix  whose  eigenvalues  have  a  reciprocal 
relationship  to  those  of  the  original  matrix.  For 
nonsingular  symmetric  matrices,  this  technique  can 
be  used  to  find  their  inverses.  We  are  looking  into 
practical  applications  of  this  development  . 
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CACHE-EFFICIENT  ALGORITHMS 
FOR  IRREGULAR  COMPUTATIONS 

Alex  Pothen 1 

Introduction  and  Motivation 

Consider  a  computation  on  an  unstructured 
mesh,  where  the  values  at  a  mesh  point  are  up¬ 
dated  using  the  current  values  at  the  neighbors 
of  that  mesh  point.  Examples  of  such  computa¬ 
tions  abound  in  scientific  computing.  In  what  or¬ 
der  should  we  visit  the  mesh  points  to  compute  the 
value  at  each  point  so  that  the  computation  is  as 
fast  as  possible? 

Observe  that  the  arithmetic  operations  (typi¬ 
cally  in  floating  point  for  scientific  computations) 
are  the  same  for  different  orders  in  which  we  visit 
the  mesh  points.  Nevertheless,  on  modern  su¬ 
perscalar  workstation  architectures,  one  ordering 
might  perform  an  order  of  magnitude  faster  than 
another  ordering.  The  reason  is  that  the  perfor¬ 
mance  of  an  algorithm  on  these  architectures  is 
limited  by  how  fast  the  memory  can  satisfy  the 
voracious  appetite  of  the  CPU  for  data.  A  fast  and 
small  memory  called  the  cache  is  used  to  feed  the 
CPU  quickly  with  the  data  that  it  needs  to  com¬ 
pute.  When  the  CPU  finds  the  data  it  needs  in  the 
cache,  then  the  computations  proceed  without  de¬ 
lay;  when  the  data  is  not  in  the  cache,  then  several 
cycles  need  to  be  spent  fetching  the  data  from  the 
primary  memory  into  the  cache.  As  modern  micro¬ 
processor  architectures  seek  to  boost  performance 
by  enhancing  the  number  of  instructions  that  can 
be  executed  in  parallel  with  pipelining,  replicated 
functional  units,  branch  prediction,  and  multiple 
issue,  it  is  crucial  to  ensure  that  the  algorithm  per¬ 
forms  the  irregular  computation  in  an  order  that 
keeps  the  number  of  cache  misses  as  small  as  pos¬ 
sible. 

Choosing  a  cache-efficient  ordering  for  an  un¬ 
structured  mesh  or  irregular  graph  is  no  simple 
matter.  In  recent  work  with  Dr.  Horst  Simon, 
Professor  Alan  George,  Dr.  Steve  Barnard,  and 
Gary  Kumfert,  we  have  described  two  ordering  al¬ 
gorithms  for  “inducing  locality”  in  irregular  meshes 
or  graphs,  and  thereby  reduce  the  number  of  cache 
misses.  One  of  these  algorithms  is  combinatorial 
in  nature,  and  is  related  to  algorithms  that  were 

1Dr.  Pothen  is  an  Associate  Professor  in  the  Computer 
Science  Department  at  Old  Dominion  University,  Norfolk, 
VA.  He  has  been  a  consultant  at  ICASE  since  1994. 


Figure  1:  A  snapshot  of  the  Sloan  algorithm  in 
progress. 

proposed  earlier  for  this  problem.  The  second  is 
an  algebraic  algorithm,  and  is  currently  the  algo¬ 
rithm  that  computes  orderings  of  the  best  quality 
for  this  problem.  Its  somewhat  intriguing  success 
can  be  justified  by  formulating  a  simplified  problem 
as  a  quadratic  assignment  problem  (a  combinato¬ 
rial  optimization  problem).  The  rest  of  this  article 
will  describe  these  results  in  more  detail. 

Wavefront  Minimization 

The  problem  considered  at  the  beginning  of  this 
article  can  be  formulated  as  a  problem  called  the 
wavefront  minimization  problem.  At  any  stage  in 
the  computation,  the  ith.  wavefront  is  the  number 
of  unvisited  mesh  points  that  are  neighbors  of  the 
i  mesh  points  that  have  been  visited.  We  seek  to 
minimize  the  sum  of  all  the  wavefronts  as  the  num¬ 
ber  of  mesh  points  visited  ranges  from  one  to  the 
total  number  of  mesh  points. 

The  wavefront  minimization  problem  has  been 
studied  in  the  sparse  matrix  algorithms  community 
because  of  its  applications  in  frontal  solvers  used  in 
structural  analysis,  and  in  incomplete  factorization 
preconditioners. 

The  Sloan  Algorithm 

A  combinatorial  algorithm  due  to  Sloan,  an  Aus¬ 
tralian  civil  engineer,  and  developed  further  by  Iain 
Duff,  John  Reid,  and  Jennifer  Scott  at  the  Ruther¬ 
ford  Appleton  Laboratory  in  Oxford,  is  currently 
the  best  in  its  class  in  terms  of  computing  the  low¬ 
est  values  of  the  total  wavefront.  This  algorithm 
takes  a  yin  and  yang  approach  to  the  wavefront 
problem,  balancing  the  conflicting  needs  of  taking 
a  global  viewpoint  of  the  mesh  while  locally  con¬ 
trolling  the  increase  in  wavefront  as  it  chooses  the 
next  meshpoint  to  be  visited. 

Sloan’s  algorithm  is  a  graph  traversal  algorithm 
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that  has  two  parts.  The  first  part  selects  a  start 
vertex  s  and  an  end  vertex  e  that  are  about  as  far 
apart  from  each  other  in  the  mesh  as  possible.  This 
is  done  by  a  few  breadth-first-searches  in  the  mesh. 
The  second  part  then  numbers  the  vertices,  begin¬ 
ning  from  s,  choosing  the  next  vertex  to  number 
from  a  set  of  eligible  vertices  by  means  of  a  prior¬ 
ity  function.  Roughly,  the  priority  of  a  vertex  has 
a  local  and  global  component:  the  local  component 
favors  a  vertex  that  increases  the  current  wavefront 
the  least,  while  the  global  part  favors  vertices  at  the 
greatest  distance  from  the  end  vertex  e. 

The  figure  shows  a  snapshot  of  the  Sloan  algo¬ 
rithm.  Seven  vertices  have  been  ordered,  and  the 
vertices  eligible  to  be  numbered  next  are  neigh¬ 
bors  of  the  ordered  vertices,  or  their  neighbors  (the 
shaded  vertices).  From  among  these  vertices,  a  ver¬ 
tex  of  largest  priority  is  chosen. 

Our  improvements  to  the  Sloan  algorithm 

Last  year,  Gary  Kumfert  and  I  provided  an 
implementation  of  this  algorithm  whose  running 
time  is  bounded  in  the  number  of  edges  in  the 
mesh;  earlier  implementations  had  time  require¬ 
ments  that  were  nonlinear  in  the  size  of  the 
problem,  0(n3/2)  for  2-dimensional  meshes,  and 
0(n5/3)  for  3-dimensional  meshes,  where  n  is  the 
number  of  meshpoints.  The  new  algorithm  is  also 
practically  faster  by  about  a  factor  of  four  on  a  set 
of  test  problems. 

We  also  varied  the  weights  of  the  local  and  global 
terms  in  the  priority  function,  and  to  our  surprise, 
found  that  our  test  problems  fell  into  two  classes 
with  different  asymptotic  behaviors  of  their  wave- 
fronts.  When  the  first  term  is  more  important  in  re¬ 
ducing  the  wavefront,  the  problem  belongs  to  Class 
1,  and  when  the  second  term  is  more  important,  it 
belongs  to  Class  2.  We  have  observed  that  the  first 
class  of  problems  represent  simpler  meshes:  e.g., 
discretization  of  the  space  surrounding  a  body,  such 
as  an  airfoil.  By  choosing  the  weights  of  the  two 
terms  in  priority  function  to  suit  each  class  of  prob¬ 
lems,  we  reduced  the  wavefront  computed  by  the 
Sloan  algorithm  by  about  a  third. 

An  algebraic  algorithm 

In  joint  work  with  Dr.  Horst  Simon  and  Dr. 
Stephen  Barnard,  I  proposed  a  novel  algebraic  al¬ 
gorithm  for  the  wavefront  minimization  problem. 
The  algorithm  associates  a  discrete  Laplacian  ma¬ 
trix  with  the  given  mesh,  and  then  computes  an 


eigenvector  corresponding  to  the  smallest  nonzero 
Laplacian  eigenvalue.  The  ordering  is  obtained  by 
sorting  the  components  of  this  eigenvector  in  in¬ 
creasing  or  decreasing  order. 

The  spectral  algorithm  can  obtain  significantly 
smaller  wavefronts  compared  to  the  combinatorial 
algorithms  for  this  problem.  The  success  of  this 
simple  algorithm  is  somewhat  intriguing.  In  re¬ 
cent  work,  Professor  Alan  George  of  Waterloo  and  I 
have  provided  stronger  justification  than  was  avail¬ 
able  for  this  algorithm. 

We  begin  by  considering  the  related  problem  of 
minimizing  the  2-sum  of  a  graph.  Given  an  order¬ 
ing  of  the  vertices  of  a  graph  (from  1  to  n.  the 
number  of  vertices),  compute  the  square  of  the  dif¬ 
ference  in  the  numbers  assigned  to  the  endpoints 
of  each  edge,  and  then  sum  this  quantity  over  all 
edges  to  obtain  the  2-sum.  The  problem  is  that  of 
computing  the  minimum  2-sum  over  all  orderings 
of  the  graph. 

We  can  formulate  the  2-sum  problem  as  a 
quadratic  assignment  problem  (QAP)  of  the  form 

min  trace  QXBXT . 

X  permutation  matrix 

where  Q  is  the  Laplacian  matrix.  B  is  the  rank-one 
matrix  with  fry  =  ij,  and  the  minimization  is  over 
the  set  of  permutation  matrices  X.  Unfortunately, 
this  problem  is  NP-complete,  and  hence  it.  is  un¬ 
likely  that  there  is  a  polynomial  time  algorithm  for 
this  problem. 

We  can  make  progress  by  minimizing  over  a 
larger  set  than  the  permutation  matrices.  We  relax 
the  requirement  that  X  have  nonnegative  elements, 
while  preserving  the  properties  of  X  being  orthogo¬ 
nal  and  having  row  and  column  sums  equal  to  one. 
This  leads  to  a  lower  bound  for  the  2-sum  problem 
by  means  of  a  projection  technique.  We  can  show 
that  the  permutation  matrix  closest,  to  the  orthogo¬ 
nal  matrix  that  attains  the  lower  bound  is  obtained 
by  permuting  the  second  Laplacian  eigenvector  in 
increasing  or  decreasing  order.  This  provides  a  rai¬ 
son  d’etre  for  the  spectral  algorithm. 

It  is  also  somewhat  surprising  that  the  lower 
bounds  obtained  are  quite  tight  (within  a  factor 
of  two)  for  the  unstructured  meshes  that  we  have 
tested.  These  lower  bounds  show  that  the  spec¬ 
tral  reordering  algorithm  can  yield  nearly  optimal 
values  of  the  2-sum,  in  spite  of  the  fact  that  mini¬ 
mizing  the  2-sum  is  an  NP-complete  problem.  To 


the  best  of  our  knowledge,  these  are  the  first  re¬ 
sults  providing  reasonable  bounds  on  the  quality  of 
the  orderings  generated  by  a  reordering  algorithm 
for  minimizing  envelope-related  parameters.  Ear¬ 
lier  work  had  not  addressed  the  issue  of  the  quality 
of  the  orderings  generated  by  the  algorithms.  We 
have  also  shown  that  problems  with  bounded  sep¬ 
arator  sizes  (n7)  have  bounded  wavefronts  (n1+7). 

Our  analysis  further  shows  that  the  spectral  or¬ 
derings  attempt  to  minimize  the  2-sum  rather  than 
the  wavefront.  Hence  a  reordering  algorithm  could 
be  used  in  a  post-processing  step  to  improve  the 
envelope  and  wavefront  parameters  from  a  spectral 
ordering.  Gary  Kumfert  and  I  have  used  the  Sloan 
algorithm  to  refine  the  ordering  produced  by  the 
spectral  algorithm  to  further  reduce  the  wavefront. 
Currently  this  algorithm  computes  the  lowest  val¬ 
ues  of  the  wavefront  on  a  collection  of  finite  element 
meshes. 

Conclusions 

Gary  Kumfert  and  I  have  designed  object- 
oriented  software  in  C++  for  the  Sloan,  RCM,  and 
the  hybrid  algorithms.  It  is  available  with  three  in¬ 
terfaces:  PETSc  (Argonne  National  Labs.),  MAT- 
LAB,  and  as  stand-alone  code. 

Current  work  involves  extending  the  Sloan  al¬ 
gorithm  for  anisotropic  problems  by  considering 
an  edge- weighted  version.  We  are  also  investigat¬ 
ing  the  application  of  another  problem,  the  1-sum 
problem,  to  the  wavefront  minimization  problem. 
We  are  also  continuing  with  more  detailed  formu¬ 
lations  of  the  cache  miss  reduction  problem. 
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